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SUMMARY 

Since turbulent duct flows play an important role in engineering, continuous efforts to investigate this 
problem have been made. Until recently, these investigations were limited mostly to experiments and various 
semi-analytic methods. Owing to the improvement of turbulence modeling, the prevailing method of prediction is 
now mainly numerical. The majority of these studies deals, however, with turbulent flows in circular or square 
ducts and only limited information is available for straight noncircular ducts. In view of this situation, we 
propose to conduct a numerical investigation of turbulent flow in a class of ducts, whose cross sections vary 
from a circle to a near square. Turbulent flow in a noncircular duct is characterized by the presence of secondary 
flow for which a more refined turbulence model than the k-e equations is required. In order to show that the 
calculated results are credible, various modes of verification were used to examine the results for a selected 
configuration, including an accuracy check by a scaling law and observing the decay of secondary flow as the 
cross section changes from a noncircular shape to a circle. After this was done, computations were performed for 
other configurations and with different Reynolds numbers from which wall shear stresses and friction factors are 
plotted. 
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contravariant velocity components (Fig. 3) 

_ . . .. — 7 — tt . _*2 — * — r ,_*2 

Reynolds stresses m cartesian coordinates, v = v /U b ,...vw = v w /U b 

cartesian coordinates, x = x*/a*...z = z*/a* 
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<{> dimensionless streamwise velocity components, <|> = U*/U* (Eq. (20)) 
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dimensional quantity 
1. INTOODUCTIOhr - 

Since turbulent flow in straight noncircular ducts is of practical importance, a wealth of investigation 
exists. Most of these studies are either experimental or semi-analytic by means of correlations through equivalent 
hydraulic diameter. See, for example, Hartnett, Koh and McComas (1962), Aranovitch (1971), and Jones (1976) 
among many others. More recently numerical computations based on turbulence models have been performed to 
predict flows in noncircular ducts. The majority of these predictions as well as experiments concern, however, 
with square or rectangular ducts. Only a relatively few are for other cross sections such as triangular, elliptical 
ducts, and rod bundles. See, for example, Carlson and Irvine (1961), Aly, Trapp, and Geirard (1978), Rapley 
(1982), and Demuren (1991). 

Measurements made by Carlson and Irvine and Aly et al. for triangular ducts demonstrate that friction 
factors for these ducts based on equivalent hydraulic diameter may be lower than for circular pipe flow. Whereas 
the use of hydraulic diameter may cause some concern, these findings still suggest that a noncircular duct may 
actually be more efficient than circular pipes. 

In many semi-analytic methods, it is a common practice to assume that there exists some form of 
similarity between a laminar and a turbulent flow for a given geometry, though the presence of the secondary 
flow is not usually considered. If we accept this assumption to be approximately correct and acknowledge the 
fact that the laminar flow friction factors for a class of ducts (super-circular ducts) are always less than those for 
pipe flow, there is a possibility that the turbulent friction factors in these ducts will be smaller than for circular 
pipe flow. On the basis of this inference, we propose to perform numerical computations for turbulent flow in 
this class of ducts. 

To predict a turbulent flow with secondary motion, the standard form of the k-e equations is no longer 
adequate and a more refined model is required. However, for the purpose of solving a fully-developed turbulent 
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duct flow, most of these models are closely related except for the model constants. See Kao (1994) for a more 
detailed discussion. Here we propose to use Speziale’s model (1987) for prediction and Naot and Rodi’s (1982) 
for comparison. 

In addition to the selection of a turbulence model, we have to generate a grid and obtain finite-difference 
equations for computation. To this end Ryskin and Leal’s method (1983) will be employed to produce an 
orthogonal and boundary-conforming grid, and the quadratic upstream differencing of Leonard (QUICK, 1979) 
will be used to solve these equations. Wall functions are also to be adopted to bridge the gap between the near- 
wall mesh and the surface. 

Because experimental data for turbulent flow in super-circular ducts are not available at present, we can 
only subject the predicted results to other forms of verification, such as the scaling law of Barenblatt (1993), the 
decay characteristics of secondary flow when a cross section changes to a circle, and the presence of multiple 
vortices without imposing symmetry conditions between every vortex. Subsequent to these verification runs, 
regular computations are carried out for other configurations, and an attempt is made to explain the observed and 
predicted fact that noncircular ducts with smooth cross-sectional contours often incur higher drag than 
corresponding ducts with sharp comers. 

2. SUPER-CIRCLE, MAPPING, AND LAMINAR FLOW COMPUTATION 


A super-circle is defined to be (see Roache 1976, p. 302) 
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where a* is the ’’radius" and n is a positive number but not necessarily an integer. For instance, when n = 2, it 
represents a circle, and as n — » ©o, it approaches to a square. See Fig. 1 for a series of configurations with 
different n’s. To construct a boundary-conforming orthogonal grid in a super circle, the weak-constraint form of 
the Ryskin and Leal method is adopted. The basic method and details are described by them and by Kao (1992), 
which transforms x and y to ^ and q. The ^ coordinate varies from zero at the origin to unity at the surface and 
is equally divided in the transformed plane. The q-coordinate corresponds approximately to the angle in a polar 
coordinate system, and varies from zero at the horizontal axis to nl 2 at the vertical axis, which is also equally 
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divided An example of such a grid is plotted in Fig. 2. This figure is for illustration; the domain for almost all 
computations is only an octant, extending from r\ = 0 to n/4. 


As mentioned earlier, the predictions of turbulence friction factors in noncircular ducts have been mostly 
empirical until recently and are mainly based on correlations through the concept of equivalent hydraulic 
diameter. Since this method lacks theoretical justification, some authors (Aranovitch and Jones, among others) 
have proposed to include laminar flow friction factor in the process of determining equivalent Reynolds number. 
Although the present investigation is numerical and does not depend on correlation, one will still find the 
availability of laminar flow solution is useful for discussion. For this reason, we make the following digression. 


The governing equation for a fully-developed incompressible laminar flow in a straight duct is 
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where all quantities are in dimensionless form and Re is a reference Reynolds number based on a reference 
velocity. See the nomenclature for its definition and Fig. 3 for other quantities. There is no secondary flow in a 
laminar flow, and as such U is the only nonvanishing velocity component and Eq. (2) is the governing equation. 
The symbol D n refers to the dimensionless pressure gradient and is a constant This is a linear equation and 
needs to be solved once for any given configuration (see Eq. (5) or (6) below). 

In order to transform Eq. (2) from x and y to £ and q, we define the following quantities: 
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in which subscriptions £, and T| refer to partial differentiations, and the metric coefficient g 12 = 0 implies 
orthogonality. Upon substitution of these quantities into Eq. (2) and after simplification, Eq. (2) becomes 
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Here U is in the axial direction and unaffected by the transformation. Notice that there is a small circle at origin 
in Fig. 2 to exclude the singularity, where Eq. (4) and mapping are invalid. This singularity is caused by the 
particular choice of £ and T| and can be removed by simply reverting it to the cartesian form locally. 

It is a relatively easy matter to solve Eq. (4) numerically, for it is basically a Poisson equation. There are 
two special cases n = 2 (circle) and n -» «> (square, Cornish 1928) for which analytic solutions are known, 
whose friction factors are 


X c = 64/Re, n = 2 
X = 56.94/Re, n -> °o 


(5) 

( 6 ) 


where Re is the actual Reynolds number based on hydraulic diameter (see the nomenclature). Similar expressions 
can be obtained numerically for other n’s. By comparing two friction factors ( XJX ) as shown in Table 1, the 
’’shape" characteristics emerge, 


TABLE 1.— COMPARISON OF FRICTION FACTORS 


n 
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2.5 

4 

6 

10 

oo 

y* 

1 

1.001 

1.012 

1.025 

1.051 

1.124 


which shows that these cross sections all have lower X’s. 

3. GOVERNING EQUATIONS AND TRANSFORMATION 

The Reynolds-averaged equations for a fully-developed incompressible turbulent flow in dimensionless 
forms are written first in a cartesian coordinate system. 
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Equation (7) is the continuity equation, Eqs. (8) to (10) are momentum equations, and Eqs. (11) and (12) the 
transport equations for k and e . To close this system, a turbulence closure model is needed. As stated earlier, 
several such models are available but they are essentially equivalent. Here we choose Speziale’s model without 
second-order terms (see Kao 1994 for justification of neglecting these terms) 
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(15) 
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There are two groups of model constants in these equations. The first group consists of the usual values 
for the standard k-e model, which are = 0.09, C el = 1.44, c e2 =i .92, o k = a e = 1.225. The second group 
C D and C E is needed for the nonlinear terms. The values chosen by Speziale are C D = C E = 1.68. In addition, 
we use Naot and Rodi’s model for calculation and comparison as mentioned previously. No equations will, 
however, be given here. 

With the metric coefficients given in Eq. (3), we transform Eqs. (7) to (18) into £-T| coordinates. In this 
process, U is regarded as a scalar, since the transformation takes place in the y-z plane. The velocity components 
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V and W are, however, to be replaced by the contravariant velocity components v and w in the £ and rj 
directions 


v = Yll W - ^V, 


w = z^V 




w 


( 19 ) 


Thus, Eqs. (9) and (10) have to be combined to form two new momentum equations for v and w. When this is 
done, a new set of equations emerges. All these equations, except those for v and w, bear close resemblance to 
their counterparts. For instance, the continuity and the U momentum equations become 
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It is not essential that Reynolds-averaged velocity components be transformed into contravariant velocity 
components. However, if V and W are retained, the system will contain both cartesian and contravariant velocity 
components, and will show up as offrdiagonal terms in the SIMPLE algorithm (Patankar and Spalding, 1972) for 
pressure correction. Whereas Shyy and Vu (1991) suggest that these off-diagonal terms may be neglected in the 
computation, we prefer to work with v and w exclusively. 

The above equations are in the high Reynolds number form. The following wall functions are needed to 
bridge the gap between the near-wall grid and the wall surface: 

U = U T In (EQ/k, w= sin \\f U T In (EQ/k 
k = U^/0.3, c = U^/Kd 

where d is the dimensionless distance in the y-z plane normal to the surface and £ is the corresponding wall 
coordinate. The symbol y is the angle between the resultant wall shear stress and the axial direction. Since this 
angle is very small, the resultant stress can be approximated by U T as in the equation for w. The above two wall 
functions for U and w are essentially Demuren and Rodi’s form (1984), including their constants E = 0.9 and 
k = 0.42, which are somewhat different from the conventional value. 
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4. NUMERICAL PROCEDURE 

The above equations after finite differencing are to be solved by means of either the QUICK method or 
Spalding’s hybrid method. The former is used for solving three momentum equations, while the latter for 
transport equations of scalar quantities. The QUICK method requires information at two grid points beyond each 
computational cell, which means that information has to be provided beyond the wall surface. To avoid this 
complication, the QUICK is implemented in the interior region two grid points away from any solid boundary, 
and the hybrid method is used to fill the gap between the interior and the wall. No change of computational 
method is, however, necessary on the other two boundaries, since these are symmetry lines. To find these 
solutions, we use the pentadiagonal matrix algorithm and a line-by-line sweep in both directions alternatively. 

Similar to laminar flow in Section 2, the transformed equations are not valid at the origin. This 
singularity can be removed by again reverting to the cartesian coordinate system locally. 

The pressure gradient dp/dx (or D n ) in Eq. (8) for a fully-developed flow can be considered as a given 
constant, but the other two gradients dp/dy and dp/dz are unknown, have to be evaluated and rectified iteratively 
by the SIMPLE algorithm. To avoid the pressure-velocity decoupling, a staggered grid is adopted for v and w, 
while permitting others to remain at the center of the control volume. 

The sequence of computation in one cycle is to solve Eqs. (9) and (10) first, the pressure-correction 
(SIMPLE) equation next, and then Eqs. (8), (11), and (12). In the iteration process under-relaxation factors are 
needed. This cycle of iteration is repeated, until the absolute value of relative residuals for difference equations 
summed over all grid points is sufficiently small, typically in the range of 1CT 5 to ltT 6 . 

Due to the limited extent of the computational domain, a relatively fine grid distribution is feasible. For 
instance, the grid distribution for Re = 6.5xl0 4 in an octant is usually 21x41 and a finer one for higher Re. One 
of the criteria to determine the grid spacing is that the distance of the near-wall mesh is £ £ 15 from the surface, 
so that wall functions can be applied properly. 

5. RESULTS AND DISCUSSION 

Unlike in the laminar flow case, the governing equations for a turbulent flow is nonlinear, and it is not 
possible to find a single solution for the entire range of the Reynolds numbers. The flow field in a given super- 
circular duct for each Reynolds number must be computed individually. Thus, a fairly large amount of 
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computations is needed to cover a moderate range of Reynolds numbers for several configurations. The 
configurations to be considered are n = 2.5, 4, and 6. 

We begin with a duct with n = 2.5. Since this configuration differs from a circle by a small amount, one 
may doubt whether a secondary flow actually exists and, if it exists, what will be its intensity. In the absence of 
experimental data, we can only provide indirect verifications. Three procedures are to be employed. First, we 
perform the computation in a half domain instead of a usual octant to see whether a similar pattern emerges. 
Secondly, the obtained result is to be subject to an accuracy check by means of a scaling law. Finally, we 
examine the decay process to see whether the secondary flow will vanish, if the cross section changes to a circle. 

5.1 Flow Field in a Half Domain 

The computational domain for regular calculations is an octant with symmetry conditions imposed along 
two "radii” (OA and OB in Fig. 3). The calculated flow field in such a domain shows clearly the presence of 
secondary flow. However, since symmetry conditions are imposed on both sides, one may still ask the question 
whether the secondary flow can subsist without these conditions. To remove this doubt, a second set of 
computation was made in which the same numerical procedure was applied except that the computational domain 
was then enlarged to include one-half of the cross section and symmetry conditions were imposed only along the 
cross-sectional symmetry line. We plot the obtained secondary flow velocity vectors after 1300 cycles of iteration 
in Fig. 4 to show that there are still four cells of clearly defined recirculatory motions. This demonstration 
suggests that the use of an octant is legitimate and a secondary flow with eight recirculatory zones is likely to be 
observed experimentally in a duct whose cross section differs from a circle by only a small amount (n = 2.5). 

5.2 Scaling law of Barenblatt 

We use this scaling law mainly to examine the accuracy of stieamwise velocity components. Although it 
resembles the conventional power-law relationship, it is formulated differendy and has the form 

= C£, a = 2 In Re/3 (20) 

where <)> = U*/U*, C is a constant, and £ is the wall coordinate. This expression is valid in the intermediate 
interval of the shear flow and implies that the graph <|> a plotted against £ is a straight line. An intermediate 
interval is defined to be a range, which is sufficiendy large in comparison with the viscous sublayer but small 
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with respect to the duct width. Strictly speaking, Eq. (20) is good only for parallel flows. Thus, using it for the 
present problem is merely an approximation with U* being only one component of the resultant. For this reason, 
the plotted graph may not be a truly straight line. However, even an approximately straight line may be 
encouraging. The coefficient a in Eq. (20) is a dimensionless empirical constant (Barenblatt and Prostokishin, 
1993). 

Plotted in Fig. 5 are four <)) a graphs for configurations n = 2.5 and 6. Only points within the intermediate 
interval are shown. It appears that groups of symbols along 0 - 15° form nearly perfect straight lines, but those 
along 0 = 0 are somewhat curved. A possible explanation is that the secondary velocity components at the center 
of vortex (0 - 15°) are nearly parallel to the surface (Fig. 4) and thus have less effect on the U* component, 
whereas those along the symmetry line are directed away from the surface, which makes the streamwise flow 
less parallel. Note also that the exponent a in Eq. (20) is a large number; a small departure from the correct 
value will cause the graph to bend. 

5.3 Decay of Secondary Motion 

A closer look at the secondary motion reveals that its intensity in a n = 2.5 duct is actually higher than 
that in a square duct at a corresponding Reynolds number. This is in contrary to our intuition. Consequently, we 
subject our computation to another test by changing the cross section from n = 2.5 to 2 (circle) and observing the 
decaying process of the secondary flow. This has to take place, if the present computational procedure is correct, 

■ n n 

owing to the fact that in a circular pipe the difference of normal Reynolds stresses v - w is zero and there is 

no mechanism to subsist a secondary motion. The process will, however, be slow, since it can only be 
annihilated through viscous dissipation. 

To illustrate this process, we plot the initial velocity profile in Fig. 6(a) and the corresponding velocity 
field based on Naot and Rodi’s model after 6000 cycles of iterations in Fig. 6(b). Whereas the remanence of the 
secondary motion is visible, most velocity vectors have essentially vanished and cannot be discerned in the 
present scale. (An open arrow head represents merely a direction, obtainable no matter how small the velocity 
vector is as long as the ratio of two components are correct, and the stem is the magnitude). As expected, the 
decay of the secondary motion is much faster in the beginning than at the later stage, and velocities slow down 
first near the surface. For instance, the flow pattern after 3000 iterations is already fairly similar to that in 
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Fig. 6(b). When the secondary flow completely subsides, the wall shear stresses in a circular pipe should be 
uniform everywhere and the friction factors agree with the Blasius or Prandtl formula. No complete subsidence is 
achieved in Fig. 6(b) especially near the plane of symmetry, and the wall shear stresses are not exactly uniform. 

As a result, the predicted friction factor is somewhat higher than the Blasius formula. This concludes our test for 
n = 2.5. 

5.4 Wall Shear Stress 

Two groups of predicted wall shear stress distributions are plotted in Fig. 7 against the dimensionless 
distance along the wall. (Wall distances are made dimensionless by referring to the duct "radius," which is a 
constant, and therefore the lengths for different configurations are different) The top group is for flows at 
relatively low Reynolds numbers, while those in the lower group are for higher Reynolds numbers. The purpose 
here is to demonstrate that cross-sectional configurations in terms of dimensionless quantities do not seem to 
have much influence on the stress distribution. As n increases, it becomes more similar to the distribution in a 
square duct In particular, it increases first, reaches the maximum at the middle of the vortex and then decreases. 
For this reason, we also plot experimental data of a square duct with similar Reynolds numbers. The resemblance 
is apparent. 

The wall stresses at one end of the vortex are always less than at the other end. This is likely due to the 
fact that the "recirculating" motion of the secondary flow at one side of the vortex (S* = 0 in Fig. 7) is directed 
away from the surface, which lifts the viscous layer and decreases the shear stress, whereas the motion at the 
other end is directed towards the surface, which compresses the viscous layer and, in turn, increases the wall 
stress. The wall shear stresses for n £ 6 increase in the first half of the vortex but bend downward in the second 
half. The reason for this is perhaps as follows. As n increases, comers begin to form, and when n is 6 or greater, 
they become relatively "sharp." A sharp comer has two walls surrounding the flow, which causes the stream wise 
velocity to slow down in its vicinity and lowers the shear stress. 

5.5 Friction Factors 

One may envision in a noncircular duct that there are two competing forces at play. One is the shape 
advantage inherent with most noncircular cross sections including various triangles (Carlson and Irvine) and the 
other is the loss induced by secondary flow showing up in the form of drag increase. When the former dominates 
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the latter as in triangular ducts, where the shape advantage can be as much as 20 percent over a circular pipe, the 
friction factors become less than those for circular pipes (see measurements taken by Carlson and Irvine and Fig. 
8). If these two forces are somehow in balance as in a square or rectangular duct, its friction factors agree with 
Blasius’ empirical equation fairly well. However, the shape advantage for a super-circular duct is relatively small 
as compared with a square duct and yet the secondary flow is still strong. Thus, the latter becomes more 
dominant and the friction factors are higher than both (see Fig. 8). 

If one considers a circle to be a curve fit for a square, which contacts the latter at four points, the 
corresponding shape for a rectangle is then an ellipse with a like aspect ratio. This geometric equivalence cannot, 
however, be translated into a friction factor equivalence, since in contrast to circular pipes the friction factors for 
elliptical ducts as measured by Cain and Duffy (1971) are considerably higher than those for rectangular duct 
(see Fig. 8). (Friction factors for rectangular ducts, though not shown in Fig. 8 to avoid cluttering, are 
approximately the same as for circular pipes) The reason is believed to be that a circular pipe, though inferior in 
shape characteristics, is absent with secondary flow, which makes it possible to have similar friction factors to a 
square duct, whereas in an elliptical duct the omnipresent secondary flow induces an additional loss and leads to 
higher friction factor than for rectangular ducts. A similar situation is believed to exist for super-circular ducts, 
which unlike a circular pipe accommodates secondary motion and suffers from higher drag than for a square 
duct. 

To extend this observation further, one may even conjecture that a duct with a smooth contoured cross 
section such as a circle, super-circle, ellipse or super-ellipse, cannot transport fluid more efficiently than a 
corresponding duct whose cross section has sharp comers. This includes both laminar and turbulent flows. Note 
that this is merely an observation and a complete justification is perhaps not attainable as our understanding 
about secondary flow losses is still inadequate (Hawthorne, 1990). 

6. CONCLUDING REMARKS 

Since there exists some evidence insofar as friction factors are concerned that noncircular ducts may be 
more efficient in transport fluid flows than circular pipes, a study was made to numerically investigate turbulent 
flow properties in ducts of super-circular cross sections. In the absence of available experimental data to test the 
validity of predicted results, we can only subject them to other forms of verification, such as an accuracy check 
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by a scaling law or observing the decaying process of the secondary flow as the cross section changes from a 
noncircular shape to a circle. The evidence seems to suggest that the predictions are correct. A somewhat 
unexpected result is that the predicted friction factors for super-circular ducts are higher than those for circular or 
square ducts. However, after a close examination of several different classes of configurations it seems that a 
noncircular duct with a smooth cross-sectional contour is likely to be less efficient than a circular pipe, but if the 
cross section is not smooth and has sharp comers, it may be as efficient as or more efficient that a circular pipe. 
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TABLE 1.— COMPARISON OF FRICTION FACTORS 


n 

2 

2.5 

4 

6 

10 

oo 


1 

1.001 

1.012 

1.025 

1.051 

1.124 
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FIGURE CAPTIONS 
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Fig. 1 Depiction of super circles of various exponents. 


Fig. 2 41x41 orthogonal grid (n = 6, equal interval in £-r| plane). 


Fig. 3 Coordinate system, velocity components and computational domain. 


Fig. 4 Vector plot of secondary motion at every other mesh point (n = 2.5, Re = 5.14xl0 4 ). 


Fig. 5 Graphs of Barenblatt’s scaling law (n = 2.5, and 6, Re = 2.18xl0 5 ). 


Fig. 6 Secondary velocity fields in a circular pipe plotted at every other point in radial direction, 
(a) Initial profile, (b) Profile after 6000 cycles of iteration at Re = 4.83xl0 4 . 


Fig. 7 Wall shear stress distributions along duct wall. 


Fig. 8 Comparison of friction factors for circular, square, triangular, elliptical and super-circular ducts. 
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Figure 1 . — Depiction of super-circles of various exponents. 
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Figure 2.-41x41 Orthogonal grid (n e 6, equal interval in plane). 
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Figure 3. — Coordinate system, velocity components and computa- 
tional domain. 
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Figure 4. — Vector plot of secondary motion at every other mesh point (n « 2.5, R© * 5.1 4x1 0 4 ). 
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Figure 5. — Graphs of Barenblatt's scaling law (n = Z5 and 6, 
R* = 2/18x105). 
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Figure 6. — Secondary velocity fields in a circular pipe plotted at every other point in radial direction, 
(a) Initial profile, (b) Profile after 6000 cycles of iteration at Rg = 4.83x1 0 4 . 
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Figure 7. — Wall shear stress distributions along duct wall. 
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Figure 8. — Comparison of friction factors for circular, square, 
triangular, elliptical and super-circular ducts. 


